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ABSTRACT 

We report an attempt to calculate energy eigenvalues of large quantum systems 
by the diagonalization of an effectively truncated Hamiltonian matrix. For this 
purpose we employ a specific way to systematically make a set of orthogonal states 
from a trial wavefunction and the Hamiltonian. In comparison with the Lanczos 
method, which is quite powerful if the size of the system is within the memory 
capacity of computers, our method requires much less memory resources at the cost 
of the extreme accuracy. 

In this paper we demonstrate that our method works well in the systems of one- 
dimensional frustrated spins up to 48 sites, of bosons on a chain up to 32 sites and 
of fermions on a ladder up to 28 sites. We will see this method enables us to study 
eigenvalues of these quantum systems within reasonable accuracy. 
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1 Introduction 

So far a lot of methods have been developed to calculate eigenvalues of large size 
quantum systems. They can be brought under two chief categories, the exact diag- 
onalization methods|I| and the Monte Carlo approaches^. In the former category, 
where the Lanczos method proves itself to be quite useful, one can expect strictly 
accurate results as long as the size of a system is small enough. The reason why 
the size of a system is limited in the Lanczos method is that it consumes enormous 
memory resources to keep every component of the state. Within this category the 
density renormalization group method || seems most prosperous to study systems 
of large sizes, but its effectivity strongly depends on both the system's properties 
and its dimensionality. By means of the Monte Carlo methods, on the other hand, 
it is possible to study quantum systems with many sites. Instead it becomes neces- 
sary to deal with the statistical errors and in some cases the so-called negative sign 
problem would prevent us to obtain meaningful results. Another disadvantage of 
the methods is that the dynamical quantities are out of their reach. 

In this paper we propose a new method to obtain the eigenvalues in quantum 
systems of large sizes, which is to be classified into the exact diagonalization ap- 
proaches. A basic idea is that we calculate expectation values 

with a choice of a trial wavefunction | and a set of trial operators {Oi} aimed to 
demand much less memory than that required in the Lanczos method, so that the 
effective Hamiltonian matrix is diagonalized within reasonable computer resources. 
A good property of our method is that it is applicable to various types of systems. 

In the next section we give a brief description of the method, giving general 
guidelines to find out the trial wavefunction and the trial operators in the systematic 
way. Concrete expressions for them will be given in the following sections together 
with the results for the spin, the boson and the fermion systems. The final section 
is devoted to summary and discussion. 
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2 Method 

In this section we outline our method to calculate the energy eigenvalues and 
their eigenstates, comparing it with the Lanczos method if necessary. 

Let us consider to calculate energy eigenvalues using a finite number of elements 
of the system's Hamiltonian H, 

Hij = (ijji | H | 
where the states {^i, fa, • • • , V'n} form a orthonormal set 

(tpi \ipj) = Sij. 

A good choice of the set is essential to obtain accurate eigenvalues using the presently 
available computer resources. 

In the Lanczos method, one begins with an orthogonal set made from an initial 

state | 

{|M/),^|M/),F 2 |vl/),...,^-i|xl/)}, (1) 

the value of n being, empirically, of order 100. Then the matrix elements are given 
by the expectation values (\& | H k \ ty) where k runs from to 2n — 1. How the 
CPU time and the memory resources needed to calculate these expectation values 
increase with the system size N? To find an answer to this question we assume, 
for concreteness, that the number of components of each state vector grows as p N 
with some positive integer p and that the Hamiltonian H is the sum of the partial 
Hamiltonians hi (i — 1, • • • , N), 

H = J2h (2) 

where hi denotes interactions between the site % and the site i + The CPU time 
to obtain the expectation values (\& | H k \ would be of order Nk provided that 
the time to calculate hi \ <fi) is in the same order for any intermediate state \ (f>). It 
should be noted that the CPU time increases not exponentially but linearly as iV 
increases. The number of the components, on the contrary, will rapidly increase as 
(pN) k / 2 if each hi produces a new state. Since one needs to keep every component of 
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the states in the evaluation, the exhaustion of the memory resources would prevent 
one to enlarge the system size N. 

Now let us introduce our method. We would like to emphasize that this method 
needs much less memory resources compared with the Lanczos method and moder- 
ate, executable CPU time to obtain results on the eigenvalues and the eigenstates 
within acceptable numerical errors. Instead of the set ([!]) used in the Lanczos method 
we employ 

{do I =| *>,6i I *),6 2 1 #>••• A-i I *)}, (3) 

where these statevectors should be linearly independent with each others and ev- 
ery operator Oi should have the same symmetry as the Hamiltonian H has — the 
translational invariance, the conservation of the total spin, the conservation of the 
particle number and so on — so that any state generated from Oi keeps the related 
quantum numbers. These requests are, however, too general to find out definite 
expressions of 0{. Being strongly motivated by the success of the Lanczos method 
we additively demand that the expectation values | H k | 1 i r ) can be, in principle, 
accurately evaluated with the set {| \l/),Oi | ^),C>2 \ ^) • • -,O n -i | Thus we 

find a hint for a systematic choice of {Oi} in the expression obtained using (0), 

N N N 

(* | H k | *> = £ £•••£<* I KK---K I *>■ (4) 

n=ii2=i ife=i 

Since the number of components one needs to calculate each | h^h^ ■ ■ ■ h ik \ ty) 
is of the order of p k ^ 2 , it seems favorable for our purpose to select the operators {Oi} 
from \hijii 2 ■ ■ ■ h ik }. Although we can deal with, of course, only a small portion of 
the whole {h^h^ • • -hi k } for systems of large sizes because the number of terms in 
the sum increases as N k , we will see it is often the case that such {Oi} is enough to 
obtain the satisfying results. 

As for the initial state | the symmetry of the state we want to study — the 
lowest energy state for example — is useful to determine it. It is desirable to choose 
a | possessing the same quantum numbers as the state to be studied. The 
restructuring method [Q for re-arrangements of states is, in addition, quite helpful 
to get the | 1 J r ) which effectively reduces the CPU time in the calculation. 
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Once the operators {Oj} and the initial state | are determined for the system 
under consideration, the rest of the calculation is conventional. We calculate the 
Hamiltonian matrix Hoij = 1, • • • , n), 

Hoij = (* | Oi-iHOj-! | (5) 

Then, employing the Gram-Schmidt method to create an orthonormal set 
{| ^i), • • • , | ip n )}, we obtain the coefficients tki, 

k 

i=i 

which are utilized to calculate the effective Hamiltonian matrix H Qi j, 

n n 

Hoij = Y Y Hoi>j>tii>tjj>, (6) 
i'=i j'=i 

where tki = for % > k. The resultant H 0i j can be easily diagonalized by means of 
the traditional numerical methods. 



3 Spin system 

In this section we apply our method to a quantum spin system. We employ 
the one-dimensional frustrated system which has both the nearest-neighbor and the 
next-nearest-neighbor interactions, concentrating our attention to its ground state. 
The reasons why we adopt this system here are that its ground state is non-trivial 
and that the exact energy eigenvalues obtained by other methods are available for 
the system of small sizes. 

The Hamiltonian we study is 

N' N' N' N' 

H — ^ &a,iVa,i+l + &b,i(?b,i+l + °a,i&b,i + <?b,i<T a ,i+l , (7) 
i=l i=l i=l i=l 

where the suffix a (b) denotes the odd (even) sites on a chain with N = 2N' spins 
and the periodic boundary condition <r a (b),N'+i = ^0(6), 1 is imposed. 
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First let us discuss on the initial state | In conventional calculations each 
state of the system is represented by | s Q) i, Sf>,i, s a ,2, $b,2, • • • , s a ,N>, Sb,N') with each 
s a n>),i = i> which indicates that the z component of the spin is ±1/2. In the 
restructuring method M we rearrange states s a ^ and into a singlet state | 0j) 
and three triplet states | 1,), | and | — 1*), 

I li) I ±0,1 > ~hi,i) i 

I — li) = I — a,i) — &,i)> (^) 

which we inclusively denote by | The complete set for the states of the system is 
now given by {| Si, S2, ■ ■ ■ , Sn>)}- Since the spin-singlet state in this set is expressed 
by only one component we employ this as the initial state, 

I *) =1 ei,e 2) ---,ejv), 

so that we can largely reduce the CPU time needed in the calculation. It should be 
minded that the contamination of the states with momentum 7r is an undesirable 
feature of this choice, because the ground state should have the zero momentum 
only. In spite of this fact, however, we will see we can obtain satisfying results with 
this I \&). 

Next we construct the operators {Oi} from H k . Taking the restructuring method 
into account, we resolve the Hamiltonian H into its components d a ,i<Jb,i and hi = 

&a,i&a,i+l + ®b.iOb,i+\ + &b,i&a,i+l, 

N' N' 
H = ^& a ,idb,i + ^2 hi. 
i=l i=l 

Then in the expansion of H k we direct our attention to the terms which are products 
of several /ij's with the suffix i being in the whole range 1 < % < N'. We combine 
those terms to express the candidates for {Oi}, which are described by 

N ' 

0(h, k 2 , ■ ■ ■ , k L ) = h i+kl h i+k , 2 ■ ■ ■ h i+kL (9) 
i=i 
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with some positive integer L, which we will call the order of the operator hereafter, 
and some positive integers k 2 , k 3 , ■ ■ -, k^. Here we keep k\ = 0. 

Now we show the numerical results. First let us report that this method presents 
the exact eigenvalues for the system of small sizes up to iV = 12. For N = 8 the set 
of the states 

{|*)Al*),..-,d B |*)} 

with 

| *) =| 01,02,03,04), 

and 

4 = 6(0), 6 2 = 6(0,1), 6 3 = 6(0, 2), 64 = 6(0,1,0), 6 5 = 6(o,i,2) 

gives us the exact eigenvalues^. 

When N = 12 the exact eigenvalue is obtainable with 28 states selected in the 
similar way. Table 1 presents the list of the operators 6, in this case. In order to 
see, without using other methods, whether the calculated eigenvalue of H is exact, 
we need to calculate the corresponding eigenvalue of H 2 , too. If the latter value 
agrees with the square of the former, the result is exact. Table 2 shows that the 
differences between the two for the N = 12 chain are negligible. 

In the cases of larger sizes we have to do with approximate eigenvalues obtained 
from a limited number of the operators. To make a systematic selection of Oi from 
its candidates, we gradually increase the value of the order of the operator L in 
Eq.(||). Suppose that we have already picked up the operators up to, say, L < L . 
The next step of the procedure is to list up 0(ki, k 2 , ■ ■ ■ , k Lo+ iYs with possible values 
of {k 2 , ks, ■ ■ ■ , k Lo+ i} [ki = 0) and add each of them into the set of {Oi} if it creates 
a new state, namely a state which is linearly independent of all the previously chosen 
Oi I ^ r )'s, when operated to the initial state | The operators in Table 1 also 
have been generated by this way. Note that the whole procedure is carried out by a 
computer program so that any work by hand is not necessary here. 

1 Note that the number of the states necessary for the set also strongly depends on the system's property. For the 
quantum spin system on a ladder where there are no cross terms {ff a i+lj f° r instance, three states are enough 
to obtain the exact eigenvalue. 
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In Table 3 we show the approximate eigenvalues of the ground state energy in 
the case of N=48 for L < 5, together with the exact one obtained by the Lanczos 
met hod 0. We see the approximation is improved as we pick up more operators 
bearing the cost of the longer CPU time, which is mainly originated from the in- 
creasing number of the states. In other words, in order to investigate systems of 
larger sizes within the available CPU time we have to adopt less precise results. So 
the maximum size of the system is determined by arranging a compromise between 
these factors. The size iV=48 is the largest one we study in this paper. 



4 Boson system 

Now let us turn our attention to boson systems with the invariant boson number 
Nfy. We consider the one- dimensional system with N sites, whose Hamiltonian is 
given by 

N i N 

H = -i ^(aj +1 a; + a\a i+1 ) + A ^(afa^) 2 , (10) 

1=1 i=l 

a,i(a\) being an annihilation (a creation) operator on the site i, with the periodic 
boundary condition a^ +1 = a[ . In the numerical work we take t = 1, while A = 4 
unless stated otherwise. 

Limiting ourselves to the case Nj, — N we employ a translational invariant state 
| 1, 1, • • - , 1), where each site of the chain has one boson on it, as the initial state 
| ^/). Note that we expect this is a good choice for large values of the coupling A 
since the state is the exact ground state in the strong coupling limit A — > oo. 

As for the operators {Oj}, we follow the same procedure as previously stated in 
Section 3 to select them from the candidates 

N 

6(k 1 ,k 2 ,---,k L ) = Y,9i+k 1 gi+k 2 ---gi+k L , (11) 

i=l 

where we use the notation cji = aja i+1 + aj +1 aj. The reason why we do not include 
the diagonal terms of the Hamiltonian in cji is that more of the inner products 
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I OiOj I \&) will become zero without the diagonal terms. 

Figs. 1 and 2 present the results on the system's ground state energy. In Fig. 1 
we plot the errors in our calculations for the N = 8 system, which is small enough 
to obtain the exact value using the Lanczos method, as a function of the coupling A. 
The results indicate that our method works well in a wide range of the parameter 
space. Especially in the strong coupling region we see the agreement is excellent, 
which confirms our expectation that the initial state we employed is good there. 
Fig. 2 shows our results on chains up to iV = 32, together with the exact values for 
iV < 10. The discrepancy between the "exact" value, which we guess from the exact 
values obtained on the short chains, and our result is about 0.5% (1%) for N = 16 
(32), respectively. These errors could be lessened by including the higher-order 
operators with L > 5. 



5 Fermion system 

This section is devoted to the study of two fermion systems, the standard Hubbard 
model and the extended Hubbard model which has been extensively studied by the 
authors using a quantum Monte Carlo method||. Here we limit ourselves only to 
the ground state of the half-filled system. We will show our method nicely works 
for these systems, too. 

The Hamiltonian of the extended Hubbard model we study is 

N r N r 

H = "^C E E J2i C h,aCi+l,l,a + h.C.} ~ t r E E[ C L<^M + h.C.} 
cr I i=l a i=\ 

N r -1 N r N r 

+ V c EE E ^M^i+l,^' +K-EE ni,a,<r n iM' + U E E n iM n i,l,h ( 12 ) 
<t,(t' I i=l a, a' i=l I i=l 

where c^} a denotes the annihilation (creation) operator for an electron, whose spin 
being a (f or |), located on the z-th rung along the leg I (a or b) of a ladder which 
has N r rungs in total and n is the number operator = c\ la Ci^ a . We assume 
the periodic boundary conditions c^ r+lla = c^\ a . The standard Hubbard model is 
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described by setting parameters V r — V c — 0. 

In the calculation we employ the restructuring method. We use a complete set 
{| a)}, where 

| a) =| Si, S 2 , ■ ■ ■ , S Nr ), (13) 

with sixteen states in Table 4 to denote the state Si on the i-th rungf]. It should be 
noted that the states in the table are eigenstates for the interaction terms between 
two sites on the i-th rung, namely eigenstates for 

o a,a' I 

The operator Oj is the product of 

h = -t c J212i c ll,* c i+hi,° + hx -} + KEE n «,i,^+i,i,^ 

a I <t,<j' I 

which comes from the interaction terms along the chains, selected from the candi- 
dates 

N 

6{ki, k 2 , ■ ■ ■ , k L ) = fi+kji+k 2 ■ ■ ■ fi+k L - (14) 

First let us report the results for the standard Hubbard model. In this case we 
choose as an initial state | *f>) a translational invariant state whose all Si are No. 9 
in Table 4. We plot the results on the ground state energy per site in Fig. 3 as a 
function of the maximum value of L, L max , for several values of N(= 2N r ) up to 
28. Here the values of the active parameters in ([12]) are set to be t c — 1, t r — 2 
and U — 4, respectively. We see that the values in Fig. 3 nicely converge into a 
value around —1.269. We would like to emphasize that this value is compatible with 
our newly calculated quantum Monte Carlo result for the same Hamiltonian with 
the different boundary conditions. Although we have to accept the poor statistical 
accuracy for the almost half-filledQ standard Hubbard model because of the negative- 
sign problem, we obtain the value —1.26 ± 0.02 at the inverse temperature /? = 10 
on an JV = 32 ladder with the open boundary conditions. 

2 Parameters u\ and 112 in Table 4 are functions of t r , V r and U given by «i = i. /l H U—V r an( j 

2 V v / (C / -^) 2 + 1 6*? 

_ 1 /, u-v r 

2 ~ 2 V V(t/-Vr) 2 + 16t?' 

3 In the Monte Carlo study the number of the fermions can slightly fluctuate according to the given value of the 
chemical potential. 
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Next we show the results on the extended Hubbard model with parameters t c — 1, 
t r = 2, V c = 2, V r = —4 and U = 4 in (|12"1). In our previous quantum Monte Carlo 
study |J it turned out that the negative-sign problem is not serious for these values 
of the parameters in a wide range of the chemical potential and that the system 
indicates a signal to the phase separation. Since we have learned from this study 
that the ground state near the half-filled system is abundant in both the hole pairs 
on the rungs and the doubly occupied sites, we choose the initial state in which the 
state on each rung % is No. 1 (No. 16) in Table 4 for even (odd) i. In Fig. 4 we plot 
the results for the ground state energy per site on the even N r ladders obtained with 
such initial states. They also show nice convergence and are close to the exact value 
—2.202 on the N = 8 ladder with the periodic boundary conditions. The quantum 
Monte Carlo result for the mostly half-filled case of the model is —2.234 ± 0.007 
(statistical error only) with (3 = 5 and N = 16 under the open boundary conditions. 



6 Summary and Discussion 

In this paper we showed an attempt to calculate eigenvalues of the Hamiltonian 
in various quantum systems of large sizes. In the calculation we start from a trial 
wavefunction, which we should choose carefully taking account of the symmetry of 
the state to be studied. Then we make a set of states by repeatedly applying some 
operators {Oi} to the trial wavefunction. Our basic idea is to determine each Oi 
systematically using several terms in the Hamiltonian of the system under consid- 
eration. Then we select a limited number of states to construct an effective matrix 
of the Hamiltonian, whose eigenvalues are easily obtained by the conventional di- 
agonalization methods. In our calculation restrictions from the memory resources 
of computers are relaxed, while the available CPU time sets bounds to the obtain- 
able accuracy and the reachable maximum system size. As concrete examples we 
demonstrate the results for low- dimensional spin, boson and fermion systems. We 
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find the method works satisfactorily for all cases we studied. 

In order to pursue further precise numerical results or to apply our method to 
larger systems including higher-dimensional ones, however, we need to achieve some 
more improvements. One of them concerns with the size of the effective Hamiltonian 
matrix. Since we select states from H k | \E'), the number of the candidate states 
rapidly increases as either the power k or the size of the system becomes large. For 
the purposes stated above, therefore, we would be obliged to construct the larger 
effective Hamiltonian matrices for which serious restrictions again would come from 
the available memory resources. So we should develop an effective method to treat 
large size matrices. Another problem is that, because of the rapidly increasing CPU 
time, it is difficult to increase the order of the operators beyond the values we have 
used here. A possible breakthrough for it would be to develop techniques to predict 
which of the states are linearly dependent. 
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Table 1: Operators 0(k\, fo, ■ ■ ■ , fcz,)'s needed for the spin system on the N = 12 chain. See (9) 
for the definition of the operator. 
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No. 


H 


H 2 


Diff. 


No. 


H 


H 2 


Diff. 


1 


-23.585 


556.24 


0.37064E-07 


15 


-7.0086 


49.121 


0.82701E-08 


2 


-22.901 


524.44 


0.37073E-07 


16 


-6.9468 


48.258 


0.25527E-07 


3 


-18.409 


338.90 


0.46285E-07 


17 


-4.2793 


18.313 


0.15361E-07 


4 


-17.271 


298.29 


0.27047E-07 


18 


-4.0000 


16.000 


0.10467E-07 


5 


-15.904 


252.93 


0.82821E-09 


19 


-3.4452 


11.869 


0.41185E-08 


6 


-13.904 


193.31 


0.29984E-07 


20 


-1.5279 


2.3344 


0.98762E-10 


7 


-13.411 


179.85 


0.44405E-07 


21 


-1.2330 


1.5204 


0.10958E-08 


8 


-12.000 


144.00 


0.26574E-07 


22 


-0.69308 


0.48036 


0.13380E-08 


9 


-11.446 


131.00 


0.42669E-07 


23 


0.47600 


0.22658 


0.14751E-08 


10 


-11.266 


126.92 


0.33501E-08 


24 


3.3595 


11.286 


0.52512E-08 


11 


-10.472 


109.67 


0.22215E-07 


25 


4.8753 


23.768 


0.89002E-08 


12 


-10.065 


101.29 


0.31402E-09 


26 


5.9206 


35.054 


0.53468E-09 


13 


-8.0485 


64.778 


0.26372E-07 


27 


10.167 


103.37 


0.16889E-07 


14 


-7.6605 


58.684 


0.26150E-07 


28 


16.676 


278.10 


0.32337E-07 



Table 2: Eigenvalues obtained in our calculations for the TV = 12 chain. Here the operators of 
the order L < 5 are used. The column of H shows the eigenvalues of the truncated Hamiltonian 
matrix, while those from the Hamiltonian squared are shown in the column of H 2 . The third value 
is the difference between the second value and the square of the first one. 
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Size 


Maximum order 


Num. of state 


E 


E( exact) 


16 


4 


37 


-1.935 


-1.94010 


20 


4 


59 


-1.928 


-1.93565 


24 


4 


92 


-1.920 


— 


28 


4 


134 


-1.911 


— 


32 


4 


191 


-1.902 


— 


36 


4 


263 


-1.892 


— 


40 


4 


352 


-1.883 




44 


4 


460 


-1.874 




48 


4 


591 


-1.865 




12 


5 


28 


-1.96538 


-1.96538 


16 


5 


78 


-1.9040 


-1.94010 



Table 3: Eigenvalues of the ground state for the truncated Hamiltonian matrix up to N < 48. 
The second column (Maximum order) indicates the maximum value of L. 
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No. 


state 


1 


1 oo> 


2 


75(4,t + c LT)l°°> 


3 


U c l, - 4t) 1 °°) 


4 


^(4a + 4i)|00) 


5 




6 




7 


^75 ( C i,a,t C lf',t + C la,J. C lb,T) 00) 


8 


7f( C !,a,T C !,a,l ~~ C Ut C U|) 00) 


9 


Nl(4a, T 4ba ~ 4aa C Lt) + M 2(4 a)T 4 aa + c\ M c\ M )\ | 00) 


10 


M^ct^ - c{ aA c\ M ) - U^cl^cl^ + 4 M 4 M )} | 00) 


11 


c IAm 1 0°) 


12 


73 ( C la,T C la,J. C lfe,T C la,T C lb,T C lfe,|) 00) 


13 


"Tf ( C !,a,t C la,J. C lb,T _ C \a,] C \,b,] C \,b,\) 00) 


14 


7f( c la,T c L,i c lfea + 4aa c !,6,T c lfea) 1 °°) 


15 


i f J J J J J J \ i nn\ 


16 


C l,a,T C l,a,J, C lb,T C lfe,J. 00) 



Table 4: States on two sites of the z-th rung used to construct a complete set. | 00) represents a 
state with no electrons on either site of the rung. Parameters u\ and ui are functions of t r , V r and 
U defined in Section 5. 
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Figure 1: Errors of the calculated eigenvalues of the ground state plotted for several values of 
the coupling A on the N = 8 chain, which are defined by the ratio (E ca i — E exact ) / E exact . The 
diamonds show the errors for the order L < 4, while the circles are for the order L < 5. 
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Figure 2: Eigenvalues of the ground state per site (the circles) calculated for TV = 8, 10, 16 and 
32 chains with the coupling A = 4 and the order L < 4. The data are plotted versus N~ 2 . The 
number of the states is 304 (1844) for N = 16 (32), respectively. Exact values for the systems up 
to N < 10 are also shown by the diamonds. 
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Figure 3: Ground state energy per site for the standard Hubbard model on ladders versus the 
maximum value of L, which is denoted by L max , calculated for several values of the system size 
N. Lines in the figure are to guide eyes. Values of the parameters in the Hamiltonian are t c = 1, 
t r = 2 and U = 4. 



19 



-2.180 




6 



Figure 4: Ground state energy per site for the extended Hubbard model described by (12) on 
ladders with size N varying up to 28, as a function of L max . Lines in the figure are to guide eyes. 
Values of the parameters in the Hamiltonian arc t c = 1, t r = 2, V c = 2, V r = — 4 and U = 4. 
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